
# Sanghyun Hong, February 23, 2019

#rm(list = ls())
###################################
## Simulation Environment Setup  ##
## Simulation Environment Setup  ##
#################################################################
#pathname<-"C:/Temp/SimulationProject/"
effectSize_List<-c(0, 50);
sigH_List<-c(0, 6.25, 12.5, 25, 50);
PubBias<-0.75;
MetaStudyN_List<- c(40, 80);
simN<-simN;
#################################################################
#
#
#
#
#
#################################################################
# Import Packages and set seed
#################################################################
library(metafor)
library(effsize)
library(sqldf) # This pakcage is for SQL
set.seed(3);
#################################################################
#
#
#
#
#
#################################################################
# Load Packages and Set Path
#################################################################
library(splines)
library(stats)
library(ggplot2)
library(gridExtra)
library(cowplot)
pathnameAK<-paste(pathname,"SimCodes/AK_Codes/",sep="")
setwd(pathnameAK)
source("Modified_AllApplications.R")
#################################################################
#
#
#
#
#
###############################################
## Primary Study Data (Cohen's d) Generation ##
## Primary Study Data (Cohen's d) Generation ##
#################################################################
PrimaryStudy_Cohen_d <- function(TreatmentEffect, sigH, primaryObs){
  obsControl<-primaryObs;
  y_cj<-rnorm(obsControl, 300, 86.603) + rnorm(obsControl, 0, 50)
  obsTreatment<-primaryObs;
  Te<- rnorm(obsTreatment, 300, 86.603) + rnorm(obsTreatment, 0, 50) + (TreatmentEffect + rnorm(1, 0, sigH))
  
  dfn1<-(length(Te)-1); dfn2<-(length(y_cj)-1);
  s1<-sum((Te-mean(Te))^2)/dfn1; s2<-sum((y_cj-mean(y_cj))^2)/dfn2;
  sp<-sqrt((dfn1*s1 + dfn2*s2)/(dfn1+dfn2))
  estimatedCohen_d<-(mean(Te)-mean(y_cj))/sp
  stdErr_Cohen_d<-sqrt(((dfn1+dfn2+2)/((dfn1+1)*(dfn1+1)))+estimatedCohen_d^2/(2*(dfn1+dfn2+2)))
  cohen_ci<-estimatedCohen_d+qt(c(0.025,.975),(dfn1+dfn2))*stdErr_Cohen_d
  cohen_t<-estimatedCohen_d/stdErr_Cohen_d;
  return(c(TreatmentEffect,sigH,TreatmentEffect/100,estimatedCohen_d,stdErr_Cohen_d,cohen_ci,(cohen_t>1.96),primaryObs))
}
###################################
## Meta Analysis Data Generation ##
## Meta Analysis Data Generation ##
#################################################################
MetaStudy <- function(effect, sigH, m, bias, PrimaryN){
  MetaStudyData<-matrix(0,nrow=m, ncol=10)
  colnames(MetaStudyData)<-c('StudyID','effect','sigH','TrueCohend','EstimatedCohend','StdErrCohend','LowerBound','UpperBound','Pos_SigEffect', 'PrimaryStudyOBS')
  if(bias==0.75){
    for(i in 1:m){
      primaryObs<-PrimaryN[i];
      if(runif(1, 0, 1)<bias){
        PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)
        while(PrimaryStudy[8]==0){PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)}
        MetaStudyData[i,]<-c(i,PrimaryStudy)
      }else{
        PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)
        while(PrimaryStudy[4]<=0){PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)}
        MetaStudyData[i,]<-c(i,PrimaryStudy)
      }}}else{
        for(i in 1:m){
          primaryObs<-PrimaryN[i]
          if(runif(1, 0, 1)<bias){
            PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)
            while(PrimaryStudy[8]==0){PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)}
            MetaStudyData[i,]<-c(i,PrimaryStudy)
          }else{
            PrimaryStudy<-PrimaryStudy_Cohen_d(effect, sigH, primaryObs)
            MetaStudyData[i,]<-c(i,PrimaryStudy)
          }}}
  return(MetaStudyData)
}
#################################################################
#################################################################
#
#
#
#
#
###################################
## Creating Tables for Output    ##
## Creating Tables for Output    ##
#################################################################
matrixRow<-length(effectSize_List)*length(sigH_List)*length(MetaStudyN_List)*simN
FE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(FE_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
RE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(RE_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP2_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP2_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
PETPEESE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(PETPEESE_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK3_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK3_Output)<-c('TrueCohend','m','sigH','PubBias','EstimatedCohend','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
SubTable<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=17))
colnames(SubTable)<-c('TrueCohend','m','sigH','FppCnt','WappCnt','WappTotal','WappSample','WappFrac','Wapp2Cnt','Wapp2Total','Wapp2Sample','Wapp2Frac','AK1P','AK3P1','AK3P2','AK3P3','I2')
#################################################################
#################################################################
#
#
#
#
#
###################################
## Simulation Begins Here        ##
## Simulation Begins Here        ##
#################################################################
start.time <- Sys.time();
sim_i<-0;
PrimaryN<-rep(c(32,64,125,250,500),100)
for(effectSize in effectSize_List){
  for(m in MetaStudyN_List){
    for(sigH in sigH_List){
      cnt<-0;
      for(simnulationLoop in 1:simN){
        cnt<-cnt+1;
        sim_i<-sim_i+1;
        data<-as.data.frame(MetaStudy(effectSize, sigH, m, PubBias, PrimaryN));
        
        ###################
        ## FIXED EFFECTS ##
        ## FIXED EFFECTS ##
        #########################################################
        effectFE<-as.matrix(data$EstimatedCohend)
        wFE<-1/data$StdErrCohend^2
        RegFE<-summary(lm(effectFE~1, weights=wFE))
        u_FE <- as.numeric(RegFE$coefficients[1])
        seFE <- as.numeric(RegFE$coefficients[2])
        p_valueFE <- as.numeric(RegFE$coefficients[4])
        ciFE <- as.numeric(confint(lm(effectFE~1, weights=wFE), level=0.95))
        covFE <- (ciFE[1]<=data$TrueCohend[1])*(ciFE[2]>=data$TrueCohend[1])
        FE_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,u_FE,seFE,u_FE/seFE,p_valueFE, ciFE[1], ciFE[2], 0, 0, 0, u_FE-data$TrueCohend[1], covFE, mean(effectFE-data$TrueCohend[1]))
        rm(effectFE, wFE, RegFE, u_FE, seFE, p_valueFE, ciFE, covFE)
        #########################################################
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## Random Effect ##
        ## Random Effect ##
        #########################################################
        effectRE<-as.matrix(data$EstimatedCohend)
        wRE<-1/data$StdErrCohend^2
        RegRE<-summary(rma.uni(data$EstimatedCohend~1,se=data$StdErrCohend, intercept=TRUE, data=data, weighted=TRUE, method="DL", level=95, digits=5));
        tau2RE <- as.numeric(RegRE$tau2)
        I2RE <- as.numeric(RegRE$I2)
        H2RE <- as.numeric(RegRE$H2)
        u_RE <- as.numeric(RegRE$b)
        seRE <- as.numeric(RegRE$se)
        p_valueRE <- as.numeric(RegRE$pval)
        ciRE <- as.numeric(c(RegRE$ci.lb, RegRE$ci.ub))
        covRE<-(ciRE[1]<=data$TrueCohend[1])*(ciRE[2]>=data$TrueCohend[1])
        RE_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,u_RE,seRE,u_RE/seRE,p_valueRE, ciRE[1], ciRE[2], tau2RE, H2RE, I2RE, u_RE-data$TrueCohend[1], covRE, mean(effectRE)-data$TrueCohend[1])
        rm(effectRE, wRE, RegRE, H2RE, tau2RE, u_RE, seRE, p_valueRE, ciRE, covRE)
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## WAAP          ##
        ## WAAP          ##
        #########################################################
        wappcnt<-1;
        u_WAAP <- as.numeric(summary(lm(data$EstimatedCohend~1, weights=1/(data$StdErrCohend^2)))$coefficients[1])
        
        APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$StdErrCohend))
        if(nrow(APSdata)<2){APSdata<-data; wappcnt<-0;}
        APSdata<-cbind(APSdata, u_WAAP=u_WAAP)
        
        wapptot <- nrow(data);
        wappsam <- nrow(APSdata);
        wappfrac <- wappsam/wapptot;
        
        effectAPS<-as.matrix(APSdata$EstimatedCohend)
        wAPS<-1/APSdata$StdErrCohend^2
        RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
        u_APS <- as.numeric(RegAPS$coefficients[1])
        seAPS <- as.numeric(RegAPS$coefficients[2])
        p_valueAPS <- as.numeric(RegAPS$coefficients[4])
        ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
        covAPS <- (ciAPS[1]<=data$TrueCohend[1])*(ciAPS[2]>=data$TrueCohend[1])
        if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-data$TrueCohend[1];}else{APS_Bias<-u_APS-data$TrueCohend[1];}
        WAAP_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
        rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## WAAP2         ##
        ## WAAP2         ##
        #########################################################
        wapp2cnt<-1;
        u_WAAP <- as.numeric(summary(lm(data$EstimatedCohend~1+data$StdErrCohend, weights=1/(data$StdErrCohend^2)))$coefficients[1])
        
        APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$StdErrCohend))
        if(nrow(APSdata)<2){APSdata<-data; wapp2cnt<-0;}
        APSdata<-cbind(APSdata, u_WAAP=u_WAAP)
        
        wapp2tot <- nrow(data);
        wapp2sam <- nrow(APSdata);
        wapp2frac <- wapp2sam/wapp2tot;
        
        effectAPS<-as.matrix(APSdata$EstimatedCohend)
        wAPS<-1/APSdata$StdErrCohend^2
        RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
        u_APS <- as.numeric(RegAPS$coefficients[1])
        seAPS <- as.numeric(RegAPS$coefficients[2])
        p_valueAPS <- as.numeric(RegAPS$coefficients[4])
        ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
        covAPS <- (ciAPS[1]<=data$TrueCohend[1])*(ciAPS[2]>=data$TrueCohend[1])
        if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-data$TrueCohend[1];}else{APS_Bias<-u_APS-data$TrueCohend[1];}
        WAAP2_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
        rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## PET-PEESE     ##
        ## PET-PEESE     ##
        #########################################################
        fppcnt <- 0
        effectFPP<-as.matrix(data$EstimatedCohend)
        wFPP<-1/data$StdErrCohend^2
        se1FPP<-data$StdErrCohend
        se2FPP<-data$StdErrCohend^2
        RegFPP<-summary(lm(effectFPP~1+se1FPP, weights=wFPP))
        ciFPP <- as.numeric(confint(lm(effectFPP~1 + se1FPP, weights=wFPP), level=0.95)[1,1:2])
        if(as.numeric(RegFPP$coefficients[1,3])>1.96){
          fppcnt <- 1
          RegFPP<-summary(lm(effectFPP~1+se2FPP, weights=wFPP))
          ciFPP <- as.numeric(confint(lm(effectFPP~1 + se2FPP, weights=wFPP), level=0.95)[1,1:2])
        }
        
        u_FPP <- as.numeric(RegFPP$coefficients[1,1])
        seFPP <- as.numeric(RegFPP$coefficients[1,2])
        p_valueFPP <- as.numeric(RegFPP$coefficients[1,4])
        covFPP <- (ciFPP[1]<=data$TrueCohend[1])*(ciFPP[2]>=data$TrueCohend[1])
        PETPEESE_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,u_FPP[1],seFPP[1],u_FPP[1]/seFPP[1],p_valueFPP[1], ciFPP[1], ciFPP[2], 0, 0, 0, u_FPP[1]-data$TrueCohend[1], covFPP,mean(effectFPP)-data$TrueCohend[1])          
        
        rm(effectFPP, wFPP,se1FPP, se2FPP, RegFPP, ciFPP, u_FPP, seFPP, p_valueFPP, covFPP)
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## A-K           ##
        ## A-K           ##
        #########################################################
        AKdata<-as.data.frame(cbind(data$EstimatedCohend, data$StdErrCohend, data$StudyID))
        attach(setup(AKdata, TRUE), warn.conflicts = FALSE)
        source("EstimatingSelection.R")
        p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
        ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
        covAK<-(ciAK[1]<=data$TrueCohend[1])*(ciAK[2]>=data$TrueCohend[1])
        
        AK_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-data$TrueCohend[1], covAK, mean(data$EstimatedCohend)-data$TrueCohend[1])
        Ak1p1<-Psihat[3]
        rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
        detach("setup(AKdata, TRUE)")
        #########################################################
        #
        #
        #
        #
        #
        #
        #
        ###################
        ## A-K 3         ##
        ## A-K 3         ##
        #########################################################
        AKdata<-as.data.frame(cbind(data$EstimatedCohend, data$StdErrCohend, data$StudyID))
        attach(setup(AKdata, FALSE), warn.conflicts = FALSE)
        source("EstimatingSelection.R")
        p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
        ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
        covAK<-(ciAK[1]<=data$TrueCohend[1])*(ciAK[2]>=data$TrueCohend[1])
        
        AK3_Output[sim_i,]<-c(data$TrueCohend[1],m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-data$TrueCohend[1], covAK, mean(data$EstimatedCohend)-data$TrueCohend[1])
        Ak3p1<-Psihat[3]; Ak3p2<-Psihat[4]; Ak3p3<-Psihat[5];
        rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
        detach("setup(AKdata, FALSE)")
        #########################################################
        SubTable[sim_i,]<-c(data$TrueCohend[1],m,sigH, fppcnt, wappcnt, wapptot, wappsam, wappfrac, wapp2cnt, wapp2tot, wapp2sam, wapp2frac, Ak1p1, Ak3p1, Ak3p2, Ak3p3, I2RE)
        rm(data)
        cat("Simulation Environment #3: ", sim_i, "/", matrixRow,"\n");
      }
}}}
#################################################################
#################################################################
#
#
#
#
#
###############################################
## Making Graphs and Saving Output           ##
## Making Graphs and Saving Output           ##
#################################################################
FE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from FE_Output group by TrueCohend, m, sigH'))
RE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from RE_Output group by TrueCohend, m, sigH'))
WAAP_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP_Output group by TrueCohend, m, sigH'))
WAAP2_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP2_Output group by TrueCohend, m, sigH'))
PETPEESE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from PETPEESE_Output group by TrueCohend, m, sigH'))
AK_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK_Output group by TrueCohend, m, sigH'))
AK3_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueCohend, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK3_Output group by TrueCohend, m, sigH'))

BIAS_FINAL <- cbind(FE_OUTPUT_FINAL$TrueCohend, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                    FE_OUTPUT_FINAL$Bias, RE_OUTPUT_FINAL$Bias,
                    WAAP_OUTPUT_FINAL$Bias, WAAP2_OUTPUT_FINAL$Bias, PETPEESE_OUTPUT_FINAL$Bias, AK_OUTPUT_FINAL$Bias, AK3_OUTPUT_FINAL$Bias)
colnames(BIAS_FINAL)<-c("TrueCohend","m","sigH","PubBias","FE_Bias","RE_Bias","WAAP_Bias","WAAP2_Bias","PETPEESE_Bias","AK_Bias","AK3_Bias")


MSE_FINAL <- cbind(FE_OUTPUT_FINAL$TrueCohend, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$MSE, RE_OUTPUT_FINAL$MSE,
                   WAAP_OUTPUT_FINAL$MSE, WAAP2_OUTPUT_FINAL$MSE, PETPEESE_OUTPUT_FINAL$MSE, AK_OUTPUT_FINAL$MSE, AK3_OUTPUT_FINAL$MSE)
colnames(MSE_FINAL)<-c("TrueCohend","m","sigH","PubBias","FE_MSE","RE_MSE","WAAP_MSE","WAAP2_MSE","PETPEESE_MSE","AK_MSE","AK3_MSE")


COV_FINAL <- cbind(FE_OUTPUT_FINAL$TrueCohend, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$coverage, RE_OUTPUT_FINAL$coverage,
                   WAAP_OUTPUT_FINAL$coverage, WAAP2_OUTPUT_FINAL$coverage, PETPEESE_OUTPUT_FINAL$coverage, AK_OUTPUT_FINAL$coverage, AK3_OUTPUT_FINAL$coverage)
colnames(COV_FINAL)<-c("TrueCohend","m","sigH","PubBias","FE_COV","RE_COV","WAAP_COV","WAAP2_COV","PETPEESE_COV","AK_COV","AK3_COV")


SubTable_Ave<-sqldf('select TrueCohend, m, sigH, avg(FppCnt) as FppCnt, avg(WappCnt) as WappCnt, avg(WappTotal) as WappTotal, avg(WappSample) as WappSample, avg(WappFrac) as WappFrac, avg(Wapp2Cnt) as Wapp2Cnt, avg(Wapp2Total) as Wapp2Total, avg(Wapp2Sample) as Wapp2Sample, avg(Wapp2Frac) as Wapp2Frac, avg(AK1P) as AK1P, avg(AK3P1) as AK3P1, avg(AK3P2) as AK3P2, avg(AK3P3) as AK3P3, avg(I2) as I2 from SubTable group by TrueCohend, m, sigH')

write.csv(BIAS_FINAL, paste(pathname,"/Output/SimulationEnvironment03_BIAS.csv", sep=""))
write.csv(MSE_FINAL,  paste(pathname,"/Output/SimulationEnvironment03_MSE.csv", sep=""))
write.csv(COV_FINAL,  paste(pathname,"/Output/SimulationEnvironment03_COV.csv", sep=""))
write.csv(SubTable_Ave, paste(pathname,"/Output/SimulationEnvironment03_Miscellaneous.csv", sep=""))


end.time <- Sys.time();
time.taken <- end.time - start.time;
cat("Total Time Taken: ", time.taken, "\n\n\n")
#################################################################
#################################################################